The Medusa of Spatial Sorting: 
Topological Construction * 
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Abstract 

We consider the simultaneous movement of finitely many 
colored points in space, calling it a spatial sorting process. 
The name suggests a purpose that drives the collection to a 
configuration of increased or decreased order. Mapping such 
a process to a subset of space-time, we use persistent homol- 
ogy measurements of the time function to characterize the 
process topologically. 

Keywords. Computational topology, persistent homology, alpha 
complexes, spatial sorting, cell segregation. 

1 Introduction 

Motivated by the observation of cell communities over time, 
we propose a topological expression of the process that facil- 
itates the identification and quantification of its features. In 
this manuscript, we focus on the mathematical framework. 

Motivation. The specific motivation for the work de- 
scribed in this paper is the experimental data on cell seg- 
regation in the developing zebrafish imaged by Heisenberg 
and Krens (9|. Making cells of two populations observable 
through fluorescent markers, they follow them through time, 
assigning each population (cell type) its own unique color. In 
this particular case, the two populations start spatially mixed 
and end in spatially segregated configurations. The segrega- 
tion is captured by a series of 3-dimensional images, which 
we turn into a shape in space-time. 
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Spatial segregation is a special case of the broader class of 
spatial sorting processes, in which we are given one or more 
distinguishable populations of particles (points in space), and 
we are interested in characterizing their spatial rearrange- 
ment in time. We aim at characterizing the spatial sorting 
process through detailed measurements of its features. The 
quantification may be used to establish a classification of 
spatial sorting processes or, on a finer scale, to differentiate 
between realizations of the same process. A common biolog- 
ical application is the establishment of phenotypes that can 
help in the classification of genetic influences. Once we have 
a description of the process beyond initial and final states, we 
may ask more subtle questions, such as whether an observed 
inverse process has a symmetric characterization. 

Results. Our contributions are primarily mathematical, 
with the goal of using the insights toward the quantitative 
analysis of experimental time-series data: 

(i) we model a spatial sorting process as a shape in space- 
time with descriptive topological properties; 

(ii) we measure this shape using the persistent homology of 
the time function; 

(iii) we provide a classification of the measurements, inter- 
preting them as aspects of the process. 

Note that measuring the process in space-time is different 
from taking the trajectory of measurements of the sequence 
of time-slices. Indeed, we will distinguish between tempo- 
rary space-time features, that can be observed in slices, from 
fleeting features that cannot be so observed. The latter re- 
quire memory and temporal reasoning and are therefore less 
readily accessible to an observer who lives in time. 

The main idea of our approach is to turn the time-series of 
geometric data into a 4-dimensional topological space whose 
connectivity is descriptive of the spatial sorting process. The 
construction takes three steps to produce the measurements: 



Step 1: construct the Voronoi tessellation to turn the data 
in a time-slice into a subset of R 3 ; 

Step 2: maintain the construction through time, effectively 
sweeping out a subset of M 3 x [0, 1]; 

Step 3 : measure the constructed subset of space-time using 
the persistent homology of its time function. 

We explain the first two steps in Section [2] and the third step 
in the Section [3] Details of the corresponding algorithms 
and their implementations can be found in 1T0| . Section [4] 
discusses the information provided by persistent homology. 
Section [5] illustrates the ideas by running the corresponding 
algorithms on simulated data. Section|6]concludes the paper. 

2 Geometry 

In this section, we introduce the sets and complexes we use 
to model spatial sorting processes. In Euclidean space, we 
call them tessellations, triangulations, and complexes, and 
in space-time, we call them medusas. They all consist of 
cells of various dimensions, but because we also talk about 
cells in the biological sense, we will use special names, such 
as Voronoi cells and simplices in Euclidean space, and stacks 
and prisms in space-time, whenever possible. 

2.1 A Moment in Time 

The input data is a finite set of colored points in K 3 . Writing 
the color as subscript, we let U be the union of U\ to Uk, 
assuming U{ n Uj = for all i ^ j. 

Voronoi tessellation and Delaunay triangulation. The 

Voronoi cell of u G U consists of all x G R 3 for which u 
minimizes the Euclidean distance among all points in U: 

vor(u) = {x G K 3 | \\x - u\\ < \\x - v\\,v G U}.(1) 

The set of Voronoi cells, V = V(U) = {vor(u) | u G U}, is 
the Voronoi tessellation of U. For each 1 < £ < k, we write 
Vi = {vor(u) | u G Ui} for the subset of Voronoi cells of 
that color. Note that V is the disjoint union of Vi , V2 , • • • , V& . 
It is often convenient to work with the dual instead of di- 
rectly with the Voronoi tessellation. Abstractly, this is the 
nerve of the set of Voronoi cells, that is: the system of sub- 
sets whose Voronoi cells have a non-empty common inter- 
section. Assuming general position of the points in U, the 
number of Voronoi cells that can have a non-empty common 
intersection is at most 4. We represent each set in the nerve 
by the convex hull of the points that generate its Voronoi 
cells, which can be a vertex, and edge, a triangle, or a tetra- 
hedron. Together, these convex hulls form a simplicial com- 
plex, known as the Delaunay triangulation of U, and denoted 




Figure 1: The full subcomplexes defined by the two colors of a 
bi-chromatic Delaunay triangulation in the Euclidean plane. 

as D = D(U). For each 1 < I < k, we write Dg for the 
full subcomplex of D that consists of all vertices of color £ 
and all edges, triangles, and tetrahedra that connect them. As 
illustrated in Figure Q~]for M 2 , the full subcomplexes of dif- 
ferent colors are disjoint, with multi-colored edges, triangles, 
and tetrahedra between them. 

Restricted Voronoi tessellation and alpha complex. In 

a loose configuration, a biological cell would generally not 
occupy the entire space alloted to it by the Voronoi cell of its 
nucleus. To better approximate the space used by the cell, 
we therefore choose a fixed positive radius, a , and consider 
the restriction of the Voronoi cell to the ball centered at the 
generating point: 

res(it) = {x G vor(u) | ||x — u|| < ao}. (2) 

Similar to before, we write R — R(U) for the set of re- 
stricted Voronoi cells, and R( C R for the subset of cells 
generated by points of color I. Clearly, R is the disjoint 
union of R±, i?2, • ■ ■ , Rk- Note that res(u) C vor(u) for 
each point u G U, It follows that the nerve of R is iso- 
morphic to a subsystem of the nerve of V. Accordingly, we 
define the dual alpha complex to consist of the simplices in 
D whose corresponding restricted Voronoi cells have a non- 
empty common intersection, denoting it as A = A(U). For 
each 1 < £ < k, we again have the full subcomplex Ag C A 
that contains all vertices of color I and all simplices connect- 
ing them. Figure|2]illustrates the definitions by restricting the 
Voronoi cells in Figure [T|to their disks. Note that At is also 
a subcomplex of D, but not necessarily a full subcomplex. 
Indeed, we have Ai = D( n A, so A# is a full subcomplex 
of D iff A e = D e . 

Homotopy equivalence. An important structural relation- 
ship between unions of Voronoi cells and their dual com- 
plexes follows from the Nerve Theorem; see e.g. ||6] page 
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Figure 2: The full subcomplexes defined by the two colors of a 
bi-chromatic alpha complex in the Euclidean plane. 



59]. To state the relationship, we ignore the difference be- 
tween a set of Voronoi cells and their union, as well as the 
difference between a simplicial complex and its underlying 
space. We write X ~ Y if the sets X and Y have the same 
homotopy type; see ifPTl page 108] for a definition. 

Lemma A. We have R e ~ At, for each color 1 < I < k. 

Proof. All cells in Re are closed and convex. It therefore 
suffices to use the corresponding special version of the Nerve 
Theorem. It says that the union of the sets and the nerve of 
the collection have the same homotopy type. The claimed 
relationships follow because Ag is a geometric realization of 
the nerves of Rg. ED 

We get Ve — De as a special case. Lemma A implies 
that Rg and Ae have isomorphic homology groups. Instead 
of dealing with the computationally more difficult union of 
Voronoi cells, we can therefore compute the homology for 
the dual complex, which is a purely combinatorial object. 



Voronoi tessellation at time t is denoted as V(t) = V(V(t)), 
and the subset of Voronoi cells of color £ is denoted as 
Vg{t) C V(t). Collecting Voronoi cells in time, we get a 
1-parameter family of cells generated by a trajectory: 

vor(u) = I^J vor(u(f)). (3) 

te[a,b] 

Noting that the Voronoi cells on the right hand side of ® 
he in distinct parallel copies of R 3 , we call vor(u) a stack. 
While the Voronoi cell in each time-slice is a 3-dimensional 
convex polyhedron, the stack itself is neither necessarily con- 
vex nor necessarily polyhedral; see Figure on the left. 
Switching fonts, we write V = V(U) for the set of stacks, 
each defined by a trajectory in U. For each 1 < I < k, we 
write Ve = Ve(U) for the subset of stacks generated by tra- 
jectories of color £. We call V the multi-chromatic and each 
Ve a mono-chromatic Voronoi medusa. 



1 


My- 




if 




P 










11 1 










f 


1 1 


i 


M 


J 


1 1 


- 





Figure 3: Left: four stack of 2-dimensional Voronoi cells. Right: 
the four corresponding prisms in the Delaunay medusa connected 
by a 3-simplex representing the flip between the two ways to trian- 
gulate four points. 



2.2 Trajectories 

Assuming a continuous trajectory for each data point, we 
form subsets of space-time by taking unions of Voronoi cells, 
both in space and in time. We begin with some definitions. 
A trajectory is a continuous mapping u : [a, b] — > R 3 , with 
< a < b < 1. We let U be a finite set of trajectories, as- 
suming no two intersect in space-time; that is: v(i) ^ u.(t) 
for all v u in U and all t G [0, 1] for which both trajecto- 
ries are defined. Writing [k] for the set {1, 2, . . . , fc}, we let 
X : U — > [k] be a coloring. At each time t € [0, 1], we have 
a finite set of points, U(i) = {u(t) | u £ U}, of course taking 
only the trajectories with a < t < b. The finite set of points 
is also colored, with coloring induced by \- 

Voronoi and Delaunay medusas. For each point u(t) £ 
U(t), we write vor(u(£)) for its Voronoi cell in R 3 x t. The 



Similar to stacks of Voronoi cells, we also consider stacks 
of Delaunay simplices, which we call prisms. There is an im- 
portant difference caused by the occasional sudden change of 
the Delaunay triangulation. Such a change is a. flip, which ei- 
ther replaces two tetrahedra by three, or three tetrahedra by 
two; see e.g. J5] page 102]. In 4-dimensional space-time, 
a flip appears as a (degenerate) 4-simplex that connects to 
the preceding Delaunay triangulations along two tetrahedra 
and to the succeeding Delaunay triangulations along three 
tetrahedra, or the other way round. Reducing the insertion 
or deletion of a point to a sequence of flips, as described 
in J5] Section 5.3], we see that prisms and 4-simplices suf- 
fice to fully describe the history of the Delaunay triangula- 
tion. We write V = £>(U) for the complex in R 3 x [0, 1] 
and T>g = T>e(\3) for the full subcomplex of color I, calling 
T> the multi-chromatic and T>e a mono-chromatic Delaunay 
medusa. Figure [3] illustrates the definitions in R 2 x [0, 1]. 
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Restricted Voronoi and alpha medusas. Following the 
distinction between unrestricted and restricted Voronoi cells, 
we extend the latter to space-time in the obvious way. We 
thus introduce a stack of restricted Voronoi cells, 

res(u) = |J res(u(t)), (4) 

t£[a,b] 

the set of such stacks, 1Z = 1Z(V), and the complex of prisms 
and 4-simplices swept out by the simplices in the alpha com- 
plex, A = A(U). Furthermore, we introduce the colored 
subsets, TZi C 7Z, and the colored subcomplexes, At C A, 
with I G [k]. We call 1Z the multi-chromatic and each TZe 
a mono-chromatic restricted Voronoi medusa. Similarly, we 
call A the multi-chromatic and each At a mono-chromatic 
alpha medusa. 

Homotopy equivalence. The structural relationship be- 
tween unions of Voronoi cells and their dual complexes ex- 
tends from M 3 to R 3 x [0, 1]. As before, we simplify the 
notation by ignoring the difference between a collection of 
cells and its union. 

Lemma B. We have Tit ~ At, for each color I e [k]. 

Proof. As in the proof of Lemma A, we appeal to the Nerve 
Theorem to give the claimed relationships. Since the cells 
are no longer convex, we use the more general version that 
applies to finite collections of closed, contractible sets whose 
intersections (of any order) are again contractible. 

While each stack in TZe of any dimension is swept out by 
a convex set and is therefore contractible, two stacks can 
intersect in two or more (lower-dimensional) stacks, which 
prevents the direct application of the Nerve Theorem. We 
finesse the difficulty by shrinking each stack at its boundary 
and thickening it into a 4-dimensional body; see Figure [4] 
There is more than one way to do this such that the bodies are 
closed and contractible, their union is the same as the union 
of the original stacks, and the common intersection of i + 1 
bodies is either empty or a contractible (4 — i) -dimensional 
set. For example, we may use the mixed complex defined 
for Voronoi cells in @ to sweep out the bodies. Letting Be 
be the resulting set of bodies, we consider the nerve, which 
we denote by Be. As illustrated in Figure |H this nerve is the 
barycentric subdivision of a complex of simplices whose ver- 
tices correspond to the original 4-dimensional stacks. Note 
that this complex is not necessarily a simplicial complex, 
since two simplices may share more than just one common 
face. Denoting this complex of simplices by At, we have 

TZe ~ Bt ~ Bt ~ At, (5) 

in which the only non-trivial, middle relation is furnished by 
the Nerve Theorem. To complete the argument, we still need 




Figure 4: Three 2-dimensional stacks (lightly shaded) giving rise to 
a decomposition into nine bodies (shaded). The nerve of the bodies 
forms the barycentric subdivision of two triangles. 

to show that At and Ae have the same homotopy type. To 
do this, we contract each prism in Ae along the time direc- 
tion. Indeed, each prism is a simplex times a time interval, 
and the contraction glues the bottom to the top face, effec- 
tively turning the prism into a simplex. We do this for all 
prisms simultaneously in order to avoid temporary inconsis- 
tencies between the prisms and their side faces, which are 
lower-dimensional prisms. Eventually, when all prisms have 
be contracted, we have turned the medusa into At. Since 
the contraction preserves the homotopy type, this implies the 
claimed relation. ED 

We get Ve — T>t as a special case. Lemma B also im- 
plies that we get isomorphic homology groups and can there- 
fore do the computation on the dual medusa, or indeed the 
complex of simplices, At, which is a purely combinatorial 
object. We refer to it as the data structure for the mono- 
chromatic alpha medusa; beyond being instrumental in the 
proof of Lemma B, it is the main computational tool we use 
to implement that algorithms that are only sketched in this 
paper; see ifTOl . 

2.3 Filtered Sequences 

We use the restricted Voronoi medusa to explain filtered se- 
quences of spaces, and free up to subscript by focusing the 
discussion on the multi-chromatic version. 

Sub- and superlevel sets. The appropriately restricted 
time function, f\n : 1Z — >• [0, 1], maps every point x E 1Z 
to its time coordinate, f\n(x). Given a moment in time, 
t € [0, 1], the corresponding sublevel set is IZt = f\^[0, t], 
and the corresponding superlevel set is 1Z 1 = /l^, 1 [t, 1]. We 
will be interested in the filtered sequences: 

C TZ C . . . C IZt C . . . C TZ Y = 1Z, (6) 
f)CK 1 C...CR t C...C^=R. (7) 
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They capture the evolution of the restricted Voronoi medusa 
in a way that allows us to compare different states without 
the burden of computing a correspondence. Similarly, we 
consider the time function restricted to the alpha medusa, 
define sub- and superlevel sets, At and A*, and form filtered 
sequences as before. Importantly, the transformation from 
the stacks to the prisms does not change the homotopy type. 

Lemma C. We have TZ t ~ At as well as TV ~ A 1 , for 
aUte [0,1]. 

We omit the proof which is similar to that of Lemma B. It 
is important to realize that Lemma C extends to all mono- 
chromatic medusas. 

Complexes of simplices. Following the proof of Lemma 
B, we find it convenient to replace the sub- and superlevel 
sets of the alpha medusa by complexes of simplices. Specif- 
ically, we write At for the complex obtained by contract- 
ing all prisms in At in time direction. Symmetrically, we 
write A 1 for the complex obtained by contracting all prisms 
in A*. There is an alternative description. Let (pbe a prism 
in A, and write f m m(<p) and / max ((^) for the minimum and 
maximum time values of the points in (p. We assign these 
two values as <7 m i n < c max to the simplex er in A obtained 
by contracting <p. Then At is the subcomplex of simplices 
with (T m in < t, and A* is the subcomplex of simplices with 
Cmax > t. Note that At and A* overlap in the simplices 
that correspond to the prisms with f m i n ((f) < t < f max ((p). 
Similar to the sub- and superlevel sets, the complexes of sim- 
plices form filtered sequences: 

C Ao C . . . C A t C . . . C Ax = A, (8) 
|ci 1 c...ci i c...ci° = i, (9) 

Again, the transformation does not affect the homotopy type. 
Indeed, the complexes forming the filtered sequences in © 
and (0 are reminiscent of the lower- and upper-star filtra- 
tions we find in (6). 

3 Algebra 

In this section, we turn the spaces of Section [2] into alge- 
braic information. The foundation of the transformation is 
the classical notion of homology, which we review. The in- 
formation is summarized in persistence diagrams, which we 
introduce for modules obtained from sub- and superlevel sets 
of the time function. 

3.1 Measuring Connectivity 

Homology groups detect and count holes in a single space. 
We begin with a brief introduction of this classical subject; 
see |[T3l for more information. 



Homology groups. Consider a simplicial complex, per- 
haps an alpha complex, A, which consists of simplices of 
dimension < p < 3. Fixing a field of coefficients, F, 
we call a formal sum of the form 7 = Ci<Ji a p-chain, in 
which the Ci are elements of F and the <7j are p-simplices in 
A, each with a fixed orientation. The boundary of the p-chain 
is the similarly weighted sum of boundaries of the simplices: 
<9 P 7 = Yl Cid p (Ti, in which d p <Ti is the sum of the (p — 1)- 
simplices that are its faces. We call 7 a p-cycle if <9 P 7 = 0, 
and we call 7 a p-boundary if there is a (p + 1) -chain 5 with 
7 = dp+18. The chains thus form chain groups connected by 
boundary homomorphisms, d p : C p — > C p _i. Similarly, the 
cycles form cycle groups, the kernels of the boundary homo- 
morphisms, Z p = kerd p , and the boundaries form bound- 
ary groups, the images of the boundary homomorphisms, 
Bp = im <9p+i. Since the boundary of a boundary is neces- 
sarily zero, we have B p C Z p and we can take the quotient, 
Hp = Zp/Bp, which is the p-th homology group. Homology 
groups can be defined quite generally, for example, as ex- 
plained above for triangulations of topological spaces. Since 
we choose the coefficients from a field, all groups we men- 
tioned are vector spaces, which are characterized by their 
ranks. For the p-th homology group, the rank is called the 
p-th Betti number, denoted as f3 p = rank H p . 

For a space of dimension k, the only possibly non-zero 
Betti numbers are fa to fa. In our case, A has dimension at 
most k = 3, and we have fa = because every 3-chain in M 3 
has non-zero boundary. The remaining three possibly non- 
zero Betti numbers have intuitive interpretations: fa counts 
components, fa counts loops, and fa counts completely sur- 
rounding walls. We get additional intuition by observing that 
the connectivity of the complement space, R 3 — A, depends 
on the connectivity of A, a relation formalized by Alexander 
Duality ifTJl p. 424]. We refer to the elements of the homol- 
ogy group of the complement space as the holes of A. Dis- 
tinguishing between the different dimensions, fa — 1 counts 
gaps between the components, fa counts tunnels passing 
through the loops, and fa counts voids surrounded by walls. 
We will compute Betti numbers for medusas in R 3 x [0, 1]. 
While more complicated than in 3-dimensional space, we 
can still interpret the Betti numbers in terms of evolutions 
of gaps, tunnels, and voids, as we will explain in Section|4] 

In addition to single spaces, we take the homology of pairs 
in which the second space is a subset of the first, such as 
(TZ,TV), for example. The p-th relative homology group, 
denoted as H p (TZ, TV ), is defined the same way as Y\ p (TV) 
except that differences in TV are ignored. In other words, two 
chains are the same if they agree on all simplices mTZ — TV, 
etc. To get an intuition for these groups, we may add the cone 
over TV. Then the rank of Y\ p (TZ,TV) is the same as the 
rank of the p-th absolute homology group for the modified 
complex, except for p = 0, where it is one less. 
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Ordinary and extended persistence modules. In the lan- 
guage of category theory, homology is a functor that maps 
a space, or a pair of spaces, to a sequence of groups, one 
for each dimension. To simplify the notation, we write H = 
(J) p Hp for the direct sum of the individual groups. This frees 
up the subscript, which we use to index the homology groups 
of a filtered sequence. Letting = t± < t2 < ■ ■ ■ < t m = 1 
be the homological critical values of fn, we note that the 
homology group is constant for all ti < t < ii+i. Hence, it 
suffices to list the finitely many groups = H (lZ ti ) to de- 
scribe the entire 1-parameter family, giving a finite sequence 
of homology groups connected by homomorphisms induced 
by inclusion: 

= R -> Ri -> . . . -> R m = H(K). (10) 

We call ([Tol l a persistence module, or more specifically the 
ordinary persistence module of /. As described in J6| Sec- 
tion VII. 1], we can use the persistence module to define 
births and deaths of homology classes and arrange them in 
pairs. Specifically, a class is born at R^ if it does not belong 
to the image of Rj_i, and it dies entering Rj if j is the small- 
est index for which the image of the class lies in the image 
of Ri_i. The birth-death pair describes a coset of homology 
classes and specifies their persistence as the absolute differ- 
ence between the function values at the two events, in this 
case tj ti. 

The ordinary module (TToT > allows for homology classes 
that are born but never die. These are the classes that de- 
scribe the entire medusa, 7Z, which are of particular interest 
to us. To get a finite measurement of duration, we extend 
the module using the relative homology groups of the form 
H(7£,7e*). Writing R m+i = H(n,n tm - i + 1 ),v/e get 

= R ^ ...-> R m -> ...-> R 2m = 0, (11) 

which we call the extended persistence module of f\n. It 
begins and ends with the trivial group, which implies that 
everything that is born also dies. We thus get a more com- 
plete set of measurements of the alpha medusa from (TTTb 
than from (ITOl . 

Isomorphic relative homology groups. Lemma C implies 
that the (absolute) homology groups of the sublevel sets of 
the restricted Voronoi medusa are isomorphic to those of the 
alpha medusa. We now show that the same is true for the 
relative homology groups. 

Lemma D. The relative homology groups H(TZ, TZ l ) and 
H(_4, A f ) are isomorphic, for all t £ [0,1]. 

Proof. Recall the exact sequences of the two pairs, which 
we write from left to right and in parallel: 

H p (A f ) ->■ H P (A) ->■ H P (A, A*)-> {A* ) -> » p -i(A) 
H p (ft*) -+H P (K)^H P (R,, H*) ^H p _! (IV) ->H p _i(ft). 



By Lemma C, the vertical maps between the first, second, 
fourth, and fifth groups are isomorphisms. The diagram 
commutes because all maps are induced by inclusion. The 
Steenrod Five Lemma thus implies that the vertical maps be- 
tween the middle groups is also an isomorphism |[T3l p. 140]. 

a 

Similar to Lemma C, Lemma D extends to the mono- 
chromatic case. 

3.2 Persistent Homology 

It is instructive to display the information contained in a per- 
sistence module as a finite multi-set of points (referred to as 
dots) in the plane. After explaining how this is done, we 
prove that the sequences of homotopy equivalent spaces in- 
troduced above give the same diagram. 

Persistence diagram. As before, we consider the re- 
stricted Voronoi medusa, 1Z, the time function f\n : 1Z —> 
[0, 1], and the extended module defined by the sub- and su- 
perlevel sets. The persistence diagram of f\n, which we de- 
note as Dgm(f\iz), is a multi-set of dots in a double covering 
of the plane; see Figure[5] Each dot has two coordinates and 
represents a coset of homology classes. These classes have a 
homological dimension, which we use to label the dot. Very 
often, we consider subdiagrams by limiting ourselves to dots 
of a particular dimension. For example, Dgm a (f\n) is the 
multiset of dots of homological dimension 0, which there- 
fore only represents components. The arrows in our draw- 
ing indicate the direction of increasing time. Since we have 
two phases, one sweeping forward and the other backward 
in time, we get every pair of coordinates twice. To explain 
this in more detail, we call the two coordinates the birth and 
the death of the dot (or of the homology classes it repre- 
sents). The birth axes appear as —45° lines in the figure, and 
the death axes as 45° lines. We decompose the persistence 
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Figure 5: From left to right: the ordinary, the horizontal, the verti- 
cal, and the relative subdiagram of the persistence diagram. The 
types of dots that can exist in the diagram of a medusa in 4- 
dimensional space-time are shown. 
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diagram into four triangles, referred to as the ordinary, hori- 
zontal, vertical, relative subdiagrams, denoted as Ord(/|-^), 
Hor(./V), Vcl(jV), Rcl(/| K ). A dot in Ord(/| K ) repre- 
sents classes that are born and die during the first phase, a 
dot in Hor(f\n) or Vc\{f\n) represents classes that are born 
during the first phase and die during the second phase, and 
a dot in Rc1(/|k) represents classes that are born and die 
during the second phase. The difference between a dot in 
Hor(/|7j) and in Vc1(/|-r.) is that for the former the birth 
coordinate is smaller than the death coordinate, while for the 
latter it is the other way round. We believe this difference is 
important, as we will explain later. 

Equivalence. Using Lemmas A to D, we can now make 
sure that the diagrams we get from the different medusas are 
the same. As before, we limit the discussion to the multi- 
chromatic restricted Voronoi and alpha medusas. Let f\n 
and f\n be the corresponding time functions. 

Lemma E. Dgm(/|_/i) and Dgia(f\-ji) are the same. 

Proof. The main tool in this proof is the diagram that con- 
nects the homology groups in the two modules: 

H(A)^-..^H(AHH(A-4 1 )^...^H(A-4°) 
H(7e )^...^H(^i)^H(^,^ 1 )^...^H(^,7e°). 

By Lemmas C and D, the vertical maps are isomorphisms. 
Since the diagram commutes, the claim follows by the Per- 
sistence Equivalence Theorem in J6j p. 159]. ED 

Similar to Lemmas C and D, Lemma E extends to the 
mono-chromatic case, and it implies the same relation for 
Delaunay and Voronoi medusas as a special case. 

3.3 Images 

We expect the measurements for the alpha medusa to be 
more meaningful than for the Delaunay medusa, but there 
is additional information in the difference. We compare by 
studying the images of the homomorphisms induced by the 
inclusion of the one in the other medusa. While we use this 
ability for mono-chromatic medusas, we simplify the nota- 
tion by describing it in the multi-chromatic case. 

Image persistence. Recall that lZ t , V}, V t , and V* are 
the sublevel and superlevel sets of the Voronoi medusas, for 
t € [0, 1]. By construction, we have TZt Q Vt and TV C V*. 
Applying the homology functor, we get finitely many dis- 
tinct groups, which we denote as Ri for the restricted and as 
Vi for the unrestricted Voronoi medusa. Aligning the two 



modules, we repeat groups as necessary and arrange them in 
a commuting diagram: 

Ro — y Ri — y . . . — y Rjv 

V -> V x -> ... -> Vjv. 

Writing : R.; — > V; for the homomorphism induced by the 
inclusion, we are interested in the sequence of images: 

= im Lq — > im ti — > . . . — > im in = 0. (12) 

Similar to the module of homology groups, (TT2T) is a se- 
quence of vector spaces connected by homomorphisms. We 
can therefore define births and deaths. We refer to the corre- 
sponding multiset of dots as the image persistence diagram, 
denoted as Dgm(im f\-ji — > /|y); see J2) for a detailed dis- 
cussion of this construction and for an algorithm. 

What does the diagram measure? In the assumed case 
of the multi-chromatic restricted included in the multi- 
chromatic unrestricted Voronoi medusa, it measures nothing 
interesting, simply because the groups are not interest- 
ing. This is different in the mono-chromatic case. Here, 
we may have a cycle defined by data points of color £ sur- 
rounding points of color different from I. In this case, we 
have a non-trivial class in a (closed or open) sublevel set 
of the restricted Voronoi medusa whose image in the cor- 
responding sublevel set of the unrestricted Voronoi medusa 
is still non-trivial. Indeed, a cycle in TZp gives rise to a dot in 
Dgm (i m f\nt ~> f | Vt ) iff it corresponds to a hole formed by 
points of color different from I. If the hole is not formed by 
such points, then it does not exist in the unrestricted Voronoi 
medusa, the class maps to 0, and there is no corresponding 
dot in the image persistence diagram. 

Instead of TZi C Vg, we can use the inclusion of IZi in 
1Z to recognize when holes are caused by interactions be- 
tween different colors. The algebraic set-up is the same, so 
we do not need to repeat it. As described in the experimen- 
tal Section [5] the latter inclusion seems to be more effective 
than the inclusion in the mono-chromatic Voronoi medusa. 
This is perhaps related to the fact that the most interesting is 
also the most difficult case in this respect, namely that of a 
1-dimensional homology class in R 3 . 

Equivalence. Similar to the conventional case, the image 
persistence diagrams do not depend on the representation of 
the space and the time function we use. Specifically, we get 
the same diagrams for the inclusion of 1Z in V as for the 
inclusion of A in T>. 

Lemma F. Dgm(im/|7 ? , — s- f\ v ) and Dgm(im /|^ -> 
f\x>) are equal as diagrams. 
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Proof. Arrange the four modules in a 3-dimensional dia- 
gram, in which the 2-dimensional section at time t consists 
of the (absolute) homology groups of the sublevel sets: 

H(A t ) -> H(TZ t ) 

I I 
H(V t ) ->■ H(V t ), 

or of the relative homology groups if the section is taken 
during the second phase. In all three directions, the maps 
are induced by inclusion, so the diagram commutes. This 
implies that we have maps between the corresponding two 
modules of images. By Lemmas C and D, the horizontal 
maps in the 3-dimensional diagram are isomorphisms, which 
implies that the maps connecting the two modules of images 
are isomorphisms. The claimed relation is implied by the 
Persistence Equivalence Theorem. ED 

Similar to Lemmas C, D, and E, Lemma F extends to the 
mono-chromatic case. We note that the two derived persis- 
tence diagrams are best computed using the complexes of 
simplices, A and T>. Indeed, these can be connected by a 
mapping cylinder, giving a complex of simplices and prisms, 
not unlike but different from the alpha medusa. From this 
complex, we get the image persistence diagrams by standard 
matrix reduction; see J2J- 

4 Classification 

In this section, we discuss the information contained in the 
extended persistence module defined by the time function, 
interpreting the corresponding events in static space-time as 
well as in dynamic temporal language. 

4.1 Plane x Time 

We begin with the 2-dimensional case as a warm-up exercise, 
but also to facilitate the comparison with the 3-dimensional 
case. Here, the medusa is embedded in M 2 x [0, 1], and the 
time function, denoted as /, maps each point to its time co- 
ordinate. We can have dots in the diagram for dimensions 0, 
1, 2, except for Ord 2 (/), Hor 2 (/), Vcl (/), and Rel (/). 
Recall that the persistence of the class represented by a dot 
X = (a, b) is the absolute difference between the two co- 
ordinates, which we denote as pers(X) = \a — b\. In the 
ordinary and relative subdiagrams, pers(X) is \[2 times the 
distance from the baseline, while in the horizontal and verti- 
cal subdiagrams, it is y/2 times the distance from the vertical 
axis that goes through the middle of the diagram. As es- 
tablished in Section [3] we have the choice between several 
filtrations, all giving the same persistence diagram. To focus 



the discussion, we limit ourselves to the medusa that best ap- 
proximates the reality of biological cells, namely the mono- 
chromatic restricted Voronoi medusa. In the text below, we 
use cell in the biological meaning of the word. 

Type analysis. The classes recorded in the ordinary sub- 
diagram are illustrated in Figure [6] (a). A dot X = (a, b) 
in Ordo(/) corresponds to a component born at t = a that 
merges with another, older component at t = b. If a = 0, 
then a non-empty subset of the cells that make up the compo- 
nent exist at the very beginning of the observed time-series. 
Its death marks the merger with another component also born 
at 0, and the tie in the decision who dies and who lives is 
broken arbitrarily. If a > 0, the component is born because 
of a new cell that is either created by cell division or en- 
ters the observation for technical reasons. A dot X = (a, b) 
in Ordi (/) corresponds to a ring (an annulus) that forms at 
t = a and whose hole fills up at t = b. 

The classes in the horizontal subdiagram are illustrated in 
Figure 0(b). A dot X = (a,b) in Horo(/) corresponds to 
a component born on the way up at t = a and dying on the 
way down at t = b. We have a < b since a is the minimum 
and b is the maximum time value of the points in the compo- 
nent. Assuming some cells of color £ live from the beginning 
to the end of the observations, we have at least one dot with 
coordinates and 1 in Horo(/). In the case of a success- 
ful segregation, in which color I eventually forms a single 
component, we have exactly one such dot and no other other 
dots with second coordinate equal to 1 in the 0-dimensional 
horizontal subdiagram. A dot X = (a, b) in Hori(/) corre- 
sponds to a ring that forms at t = a and breaks up at t = b. 
The breaking up of the ring is not detected by the homology 
of the growing sublevel set, but rather later by the relative 
homology during the down phase. The classes recorded in 
the vertical subdiagram are illustrated in Figure|6](c). A dot 
X = (b, a) in Vcli(/) corresponds to a tunnel born on the 
way up at t = b which dies on the way down at t = a, where 
a < b. The tunnel is the plane-time expression of a gap that 
temporarily opens up between two subpopulations of a com- 
ponent that later closes again. A dot X = (b, a) in Vci2(/) 
corresponds to a void born on the way up at t = b which 
dies on the way down at t = a, where a < b. The void is 
formed by a temporary ring structure we see within an inter- 
val of level sets. In contrast to the rings discussed earlier, this 
ring is formed by puncturing. First, the hole expands until it 
eventually shrinks back to a point and disappears at t = b. 

The classes recorded in the relative subdiagram are illus- 
trated in Figure |6](d). A dot X = (b, a) in Rcli(/) corre- 
sponds to a component of the level set that splits off another 
component at t = a, and dies off at t = b. Since this event 
is detected in the down phase of the module, we see the birth 
of a 1-dimensional class at t = b and its death at t = a. A 




(a) Ordinary subdiagram. Left: a gap appears when a new component (b) Horizontal subdiagram. Left: a gap appear and disappears when a new 

is born and disappears when this component merges with another com- components is formed and ends. Right: a ring appears when two ends of 

ponent. Right: a ring appears when two ends of a component meet and a component meet and disappears when the ring breaks open, 
disappears when the surrounded hole fills up. 




(c) Vertical subdiagram. Left: a gap appears when a component splits (d) Relative subdiagram. Left: a gap appears when a new component is 

in two and disappears when the two merge back together into a single split off and disappears when the component ends. Right: a ring appears 

component. Right: a ring appears when a disk is punctured, which creates when the disk is punctured and disappears when the ring breaks open, 
a hole that first grows and later fills up. 



Figure 6: The classes recorded in the four subdiagrams of the persistence diagram. 



dot X = (b,a) in Rci2(/) corresponds to ring formed by 
puncturing a disk at t = a. The hole expands and eventually 
disappeared because of the breaking up of the ring at t = b. 

Interaction between colors. The appearance and disap- 
pearance of a hole may or may not be facilitated by in- 
teractions with cells of different color. For example, a 
gap recorded in the ordinary subdiagram may disappear by 
squeezing out cells of color different from £, or simply by 
locally consolidating the configuration. In the former case, 
it is likely that the gap also existed in Vi and unlikely that 
it existed in 1Z. Accordingly, the gap in IZe would be repre- 
sented by a dot in the image persistence diagram of IZe Q Ve 
but not in that of IZg C 1Z. In the latter case, it would of 
course be the other way round. Similar distinctions can be 
made for gaps recorded in Horo, Vcli, and Rcli. The analy- 
sis of rings formed by cells of color I is similar. A ring may 
form around a population of cells of color different from I, 
or around empty space. In the former case, the ring can only 
disappear by killing the surrounded cells or by breaking up 



to release these cells. Since killing seems unlikely, this may 
mean that such rings are rare in the ordinary and more com- 
mon in the horizontal subdiagram. Similarly, the latter case 
is likely to be the preferred configuration for rings recorded 
in the vertical and the relative subdiagrams. Indeed, these 
rings form by puncturing and can therefore not be invaded 
by cells from the outside. Note, however, that the suggestion 
that a dot in the persistence diagram of IZe either appears or 
does not appear in the image persistence diagram is an over- 
simplification. Thinking of a dot as an interval, it can break 
up into shorter intervals, or merge with others into a longer 
interval. Similarly, there is generally no simple relationship 
between the image persistence diagram and the diagram of 
the containing medusa. 

4.2 Space x Time 

Similar to the planar case, we discuss the time function on a 
mono-chromatic restricted Voronoi medusa, [0, 1] 

with 7^ C R 3 x [0, 1], 
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Types. Compared to the planar case, we have one more 
dimension and therefore one additional case in each subdi- 
agram; see Figure [5] In temporal language, we follow the 
evolution of 0-, 1-, and 2-dimensional holes, which we refer 
to as gaps, tunnels, and voids. There are two ways a hole 
can appear: by closing the surrounding cycle, or by punc- 
turing. Of course, a hole can already exist at the beginning, 
at t = 0.0, in which case its birth coordinate is 0.0. Simi- 
larly, there are two ways a hole can disappear: by breaking 
up the surrounding cycle, or by filling up the hole. In ad- 
dition, it may remain to the end, at t = 1.0, in which case 
the death coordinate is 1.0. We observe that a p-dimensional 
hole in the level set gives rise to either a p-dimensional or 
a (p + l)-dimensional class in the persistence module, and 
which it is depends on how the hole comes into existence: 
either by closing a p-dimensional cycle or by puncturing a 
(p + 1) -dimensional chain; see Table Q] A hole recorded 
in the ordinary subdiagram is formed by closing a cycle, al- 
lowing for the case that it exists already at the beginning, 
and it disappears by filling up. We refer to this action as 
aggregation. Symmetrically, a hole recorded in the relative 
subdiagram is formed by puncturing, and the surrounding 
cycle either breaks up or the hole remains until the end. We 
refer to this action as disaggregation. In contrast, the holes 
recorded in the vertical and horizontal subdiagrams represent 
aggregations and disaggregations that are either transient or 
incomplete. Specifically, a hole recorded in the vertical sub- 
diagram appears by puncturing, but instead of breaking up, 
it disappears by filling. Finally, a hole in the horizontal sub- 
diagram forms by closing a cycle that breaks up later, which 
allows for the case that the hole exists already at the begin- 
ning or that is remains to the end. 





Ord 


Hor 


Vcl 


Rel 




close/fill 


close/break 


punct/fill 


punct/break 





gap 


gap 






1 


tunnel 


tunnel 


gap 


gap 


2 


void 


void 


tunnel 


tunnel 


3 






void 


void 



Table 1: The twelve cases that arise from four kinds of evolutions 
possible for each of three different types of holes in a 1-parameter 
family of spaces in R 3 . 



Rings and tunnels. The 0- and 2-dimensional holes in R 3 
behave like the 0- and 1-dimensional holes in M 2 , while the 
tunnels introduce behavior that cannot be observed in the 
planar case. Drawing the ring passing around a tunnel as 
a solid torus, Figure [7] shows the four different evolutions 
recorded by dots in the four subdiagrams of Dgm(/). 

As in the other cases, it is interesting to distinguish be- 
tween tunnels formed by interactions with cells of color dif- 




Ordj Hor-L Vcl 2 Rel 2 



Figure 7: From left to right: the evolution of a ring/tunnel with 
recording dot in the ordinary, horizontal, vertical, relative subdi- 
agram of Dgm(/). The critical configurations in the second and 
fourth rows separate the regular configurations in the other rows. 

ferent from I, from tunnels formed without such interactions. 
As before, we use the image persistence diagrams induced by 
TZe C Vi or Hi C TZ. If the tunnel corresponds to a dot in the 
diagram induced by the inclusion in Ve, then this is evidence 
for a different color cell squeezing through an opening in the 
wall of color I. Similar evidence for this event is provided 
by the absence of this dot in the diagram induced by the in- 
clusion in TZ. The different color cell may be successful and 
move from one side of the wall to the other, or it may be 
repelled by the force holding the wall together and bounce 
back. Which case it is may sometimes be apparent but can 
generally not be read off the diagrams. 

5 Experiments 

This section presents results for datasets generated with a 
software that simulates the spatial motion of cells. It serves 
as an illustration of our topological methods, and as a proof 
of concept. The results encourage us to proceed to the orig- 
inal goal of applying the method to cells of developing ze- 
brafish embryos in the near future. 

Cell segregation. Cell sorting involves both the segrega- 
tion of a mixed population of cells with different properties 
into distinct domains, and the active maintenance of their 
segregated state. It has been described to occur in vivo in 
a wide variety of biological processes, such as the early em- 
bryonic development of multi-cellular organisms. The segre- 
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gation behavior can be studied in vitro - in artificial cell cul- 
ture conditions - by mixing cells with different properties, 
and observing their subsequent sorting behavior. Such stud- 
ies have been instrumental in revealing the cellular properties 
driving the segregation process. Different theories have been 
postulated and experimentally tested to explain the mecha- 
nism that drive cell sorting both in vitro and in vivo. How- 
ever, to determine relative contributions of cellular properties 
requires the accurate recording as well as a quantitative and 
descriptive analysis of the sorting process in space and time. 
This is the purpose of this paper, namely a refined measure- 
ment of the spatial process that will allow us to compare dif- 
ferent experimental conditions and quantify similarities and 
differences to wild type processes. 

Here we use the zebrafish gastrulation as a case-study, in 
which mesoderm cells are induced and segregate from the 
overlying ectoderm germ layer. This process can be recapit- 
ulated in vitro by randomly mixing population of ectoderm 
and mesoderm cells in culture, and using advanced time- 
lapse microscopy imaging techniques ifTTl to observe how 
mesoderm cells engulf the ectoderm cells. 

Simulations. Because of the difficulties of accurately 
tracking the trajectories of cell nuclei in real data, we concen- 
trate on simulated data to demonstrate our methods. We use 
the publicly available simulation software CompuCelUI^ 
lfl4l to get data sets that imitate the cell segregation process. 
The software uses a 3-dimensional version of the widely 
applied Cellular Potts model (also known as Glazier and 
Graner model JS]) to describe in vitro segregation of ecto- 
derm and mesoderm cells. Each cell is represented by a set 
of voxels (unit integer cubes in E 3 ), and randomly tries to 
extend into the empty space or a neighboring cell. An ele- 
mentary change is the alteration of the membership status of 
a single voxel. The decision whether or not to accept such a 
change depends on an energy function that encodes the char- 
acteristics of the dynamic process. Evaluating the energy 
before and after the change, the difference is translated into 
the probability of acceptance. It is higher for negative differ- 
ences, which drive the configuration toward smaller energy. 

We give more details on our specialization and refer to 
lfl4l for the general simulation framework. The first term of 
our energy function constrains the shape of the cells. Setting 
the target volume to 49.0, each cell increases the energy by 
the square of its deviation from the target value. The same 
holds for the surface area, where we set the target to 26.0. 
The software allows for weights to control the relative influ- 
ence of these terms, and we use weights 100.0 for the volume 
and 3.0 for the area. The second term of our energy function 
influences the inter-cellular behavior. We call two neighbor- 
ing voxels an interface if they belongs to different cells, or 

1 http://www.compucell3d.org/ 



to a cell and empty space. In our setup, an interface between 
two cells decreases the energy, so that cells tend to stick to 
each other. Specifically, we define two cell types, which we 
call red (simulating ectoderm) and blue (simulating meso- 
derm). An interface between two blue cells, or between a red 
and a blue cell decreases the energy by 200.0, while an in- 
terface between two red cells decreases the energy by 400.0. 
This results in red cells sticking to each other stronger than 
blue cells. Finally, an interface between a cell and empty 
space increases the energy by 20.0 for a blue cell, and by 
100.0 for a red cell, with the effect that cells generally avoid 
contact with the surrounding space, and red cells do so with 
more determination than blue cells. 

Datasets. We construct the initial configuration of our ex- 
periments by tiling the domain into 6x6x6 = 216 cubes 
and placing one cell into each cube to fill it completely. This 
setup leads to cells that are larger than the target volume, re- 
sulting in a high initial value of the energy function. In the 
first steps of the simulation, the cells shrink and cluster in 
the middle of the available space. However, some of them 
lose contact and stay separated from the bulk for a long time, 
if not the entire time. We consider three simulations. In the 




Figure 8: The restricted Voronoi tessellation at four moments in 
time. At the beginning, the cells form a cubical grid (upper-left). 
Next, the cells move toward the middle of the available space, leav- 
ing some outliers behind (upper-right). Thereafter, the blue cells 
begin to engulf the red cells (lower-left). Finally, the blue cells 
form a sphere surrounding a ball of red cells (lower- right). 

first, we color each cell either red or blue according to a fair 
coin flip; see Figure[8] This intermixed configuration mimics 
the biologically interesting process of cell segregation. The 
second and third are control experiments in which all cells 
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are of the same type, either red or blue. To obtain trajecto- 
ries, we compute the center of gravity of each cell at various 
moments in time and interpolate them piecewise linearly. In 
all experiments, we set the radius of the ball approximating 
one cell to ccq = 4.0, and we model the arrangement of cells 
by the thus defined restricted Voronoi tessellation. 

As explained in the preceding sections, the behavior of the 
cell population in time is characterized by the persistence and 
image persistence diagrams of the time function restricted to 
the various medusas. There are too many diagrams to fit in 
this paper, so we show a subset, selected to illustrate inter- 
esting behavior of the data. For each diagram, we count the 
dots and we compute the 1-norm, which we define as 

UDgmC/Ollj = Yl P ers ( X )' < 13 > 

XeDgm(f) 

These numbers are collected for each dimension and each 
subdiagram separately, giving a good overall impression of 
the activity recorded in the diagram. 




Figure 9: The 2-dimensional persistence diagrams of the restricted 
Voronoi medusa of B. 



Multi-chromatic data. We start with the discussion of the 
intermixed simulation, letting R and B denote the corre- 
sponding sets of red and blue cells. See Figure [8] for il- 
lustrations of the restricted Voronoi tessellation, Figure [9] 
for one of the persistence diagrams, and Table [2] for the 
quantitative summaries of all diagrams, divided up into di- 
mensions and subdiagrams. The only high-persistence fea- 
ture in R is recorded in the horizontal subdiagram and sim- 
ply represents the component itself. There are additional 
high-persistence features recorded in the 1-dimensional ver- 
tical and relative subdiagrams of B, representing the isolated 
(blue) cells that lose contact to the bulk during the initial 
phase. Cells recorded in the vertical subdiagram merge back 
into the bulk, while cells recorded in the relative subdiagram 
remain separate until the end of the simulation. In addition, 
the 2-dimensional horizontal subdiagram of B records fea- 
tures whose persistence adds up to a value of 0.76. They 
describe the sphere of blue cells that wrap around the red 



cells. The fact that it is represented by several instead of a 
single dot indicates that the sphere is pinched from time to 
time during the simulation; compare with Figure [8] 

To cast additional light on the process, we show results 
for the blue restricted Voronoi medusa included in the blue 
unrestricted Voronoi medusa; compare the two bottom pan- 
els of Table [2] We no longer see any persistent features in 
the vertical or relative subdiagrams because the outliers are 
merged with the bulk by the unrestricted medusa. Further- 
more, the image persistence diagram features a single hori- 
zontal void with high persistence, which is the consolidation 
of the many low-persistence horizontal voids of the blue re- 
stricted Voronoi medusa. 



R 


Ord 


Hor 


Vcl 


Rel 


sum 


gaps 
tunnels 
voids 


17 0.03 
103 0.04 
21 0.00 


1 1.00 
5 0.00 
0.00 


7 0.03 
0.00 
2 0.00 


0.00 

1 0.00 
0.00 


25 1.07 
109 0.05 
23 0.00 


sum 


141 0.08 


6 1.00 


9 0.03 


1 0.00 


157 1.12 


B 












gaps 
tunnels 
voids 


10 0.02 
237 0.27 
32 0.00 


1 1.00 
37 0.05 
61 0.76 


22 1.58 
21 0.04 

0.00 


3 2.55 
41 0.02 

0.00 


36 5.15 
336 0.40 
93 0.76 


sum 


279 0.30 


99 1.82 


43 1.62 


44 2.57 


465 6.32 


B, Vor 












gaps 
tunnels 
voids 


0.00 
9 0.02 

1 0.00 


1 1.00 
4 0.01 
1 0.95 


0.00 
3 0.00 
0.00 


0.00 
7 0.00 
0.00 


1 1.00 
23 0.04 

2 0.95 


sum 


10 0.02 


6 1.96 


3 0.00 


7 0.00 


26 2.00 



Table 2: From top to bottom: the sizes and 1-norms of the persis- 
tence diagrams of the red and blue restricted Voronoi medusas, and 
of the image persistence diagrams defined by the inclusion of the 
blue restricted Voronoi medusa in the blue Voronoi medusa. Bold- 
face numbers correspond to dots in Figure[9] 




Figure 10: The 1-dimensional persistence diagrams of the restricted 
Voronoi medusa of BB. 



Mono-chromatic data. The control experiments use sim- 
ulations with cells of a single type, which is either red or 
blue. To avoid confusion with the data from the previous ex- 
periment, we denote the mono-chromatic sets of cells by RR 
and BB. Table [3]presents the quantitative summaries of the 
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topological analysis. The process looks similar for both col- 
ors: there is a component formed by the majority of the cells 
(represented by a dot in horizontal subdiagram), and there 
are several outlier cells (represented by dots in the vertical 
and relative subdiagram). Moreover, many low -persistence 
ordinary tunnels are created in the initial phase, whose per- 
sistence adds up to a significant number. 

There are also differences. First, the blue cells create 
many vertical gaps (components that split and re-merge dur- 
ing the simulation), while they are rare for the red cells. An 
inspection of Figure [K)]reveals that most of these splits hap- 
pen early in the simulation, when the cells still shrink. This 
effect can be explained by the smaller stickiness of the blue 
cells, which makes it more likely for them to lose contact 
with the bulk. Once separated, outliers move randomly in 
space and find contact with the bulk only by chance. Curi- 
ously, the red medusa has the same number of relative gaps 
(components that split and the pieces stay separate) as the 
blue medusa. It suggests that the red outliers find it more 
difficult to merge back into the bulk. One reason may be 
their stickiness to each other, and their phobia of the empty 
space, which limits their search activity. 

A second difference can be observed in the vertical voids, 
which only exist for the red cells. Recall that they are created 
by puncturing and destroyed by filling. As we can read off 
the table, these holes are short-lived. There is no biological 
reason why cells would create such holes, and they should 
indeed be considered an artifact of the abstraction of cells as 
relatively inflexible restricted Voronoi regions. Nevertheless, 
even this artifact tells us something about the simulation pro- 
cess. The presence of such holes means that the red cells tend 
to deviate more form the shape of a round ball than the blue 
cells. This is consistent with the energy function in our setup 
in which red-red interfaces are more favorable than blue-blue 
interfaces. Therefore, the red cells get more reward for in- 
creasing their interface area at the expense of deviating from 
the target volume and surface area. 



RR 


Ord 


Hor 


Vcl 


Rel 


sum 


gaps 
tunnels 
voids 


0.00 
607 0.60 
339 0.03 


1 1.00 
47 0.04 
0.00 


2 0.01 
1 0.00 
24 0.02 


3 2.99 
3 0.00 
10 0.00 


6 4.01 
658 0.64 
373 0.06 


sum 


946 0.64 


48 1.04 


27 0.04 


16 2.99 


1037 4.72 


BB 








gaps 
tunnels 
voids 


0.00 
582 0.55 

310 0.02 


1 1.00 
81 0.09 

0.00 


15 1.52 

0.00 
0.00 


3 2.56 

16 0.00 

4 0.00 


19 5.08 
679 0.65 
314 0.02 


sum 


892 0.58 


82 1.09 


15 1.52 


23 2.56 


1012 5.76 



Table 3: The sizes and 1-norms of the persistence diagrams of the 
red medusa at the top and the blue medusa at the bottom. Boldface 
numbers correspond to dots in Figure [Tol 




(a) Image diagram induced by the inclusion of the restricted Voronoi 
medusa of B i in unrestricted Voronoi medusa of the same set of cells. 




(b) Image diagram induced by the inclusion of the restricted Voronoi 
medusa of Bi in the restricted Voronoi medusa of BB. 



Figure 11: The 1-dimensional image persistence diagrams of the 
inclusion of Bi in its corresponding Voronoi medusa at the top, 
and the blue combined medusa at the bottom. 

Image diagrams. We create a randomly intermixed pro- 
cess by decomposing BB into two sets of cells, B x and B 2 , 
with fair coin flips. There is no reason to expect a sorting 
process similar to the simulation with red and blue cells. In- 
stead, the two (identical) cell types form a random pattern 
with arbitrary neighbor exchanges over time. 

The quantitative results for Bi are given in Table HJ those 
for B2 are similar. Note the large number of tunnels, which 
are distributed over the entire simulation. Recall that we 
have two possibilities for computing image persistence: by 
mapping the restricted Voronoi medusa of Bi into the un- 
restricted Voronoi medusa of the same set of cells, or into 
the restricted Voronoi medusa of BB. For comparison, we 
display the quantitative results of both inclusions in Table @] 
showing the corresponding 1-dimensional diagrams in Fig- 
ure QT| The majority of features - both in terms of num- 
ber and persistence - carry over to the unrestricted Voronoi 
medusa, which suggests that most of the holes in the medusa 
of Bi are caused by the presence of cells in B 2 . This is 
consistent with the relative sparsity of the image diagrams 
defined by the inclusion of Bi in BB. Notable exceptions 
are the gaps recorded in the 1-dimensional vertical and rela- 
tive subdiagrams, which describe the outliers separated from 
the bulk mentioned earlier. 
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Bi 


Ord 


Hor 


Vcl 


Rel 


sum 


gaps 
tunnels 
voids 


1 1 0.00 
1406 4.32 
32 0.00 


1 1.00 
629 5.41 
0.00 


42 1.28 
648 6.02 
0.00 


1 0.76 
1273 4.60 
0.00 


55 3.05 
3956 20.36 
32 0.00 


sum 


1449 4.33 


630 6.41 


690 7.31 


1274 5.36 


4043 23.42 


Bi, Vor 












gaps 
tunnels 
voids 


0.00 
1113 4.0S 

0.00 


1 1.00 
445 3.70 

2 0.00 


10 0.66 

627 5.92 
0.00 


0.00 

1183 4.47 
0.00 


11 1.66 
3368 18.15 
2 0.00 


sum 


1113 4.05 


448 4.70 


637 6.58 


1183 4.47 


3381 19.82 


Bi, BB 












gaps 
tunnels 
voids 


0.00 
149 0.11 

32 0.00 


1 1.00 
15 0.01 

0.00 


15 1.04 

0.00 
0.00 


1 0.76 

1 0.00 
0.00 


17 2.80 
165 0.13 
32 0.00 


sum 


181 0.11 


16 1.01 


15 1.04 


2 0.76 


214 2.94 



Table 4: From top to bottom: the sizes and 1-norms of the per- 
sistence and image persistence diagrams of the first blue medusa. 
Boldface numbers correspond to dots in FigurefTTl 



6 Discussion 

This paper introduces the persistent homology analysis of a 
medusa as a novel method to measure cell segregation. To 
provide a proof of concept, we have computed these mea- 
surements for simulated time series of 3D data. The appli- 
cation of the method to imaged cell segregation in zebrafish 
embryos is forthcoming. 

The medusa introduced in this paper is related to the vine- 
yard described in |6] Section VIII. 1]; see also (3). How- 
ever, there are differences. Specifically, the vineyard would 
be constructed for two parameters, the restricting radius, a, 
and the time, t. The result is a richer structure, namely a 
collection of curves in R 2 x [0, 1], which is therefore more 
difficult to understand. In contrast, the medusas fix a to ao 
and in this way facilitate a more compact representation of 
a subset of the vineyard. This is appropriate for biological 
cells whose size does not vary substantially with time, and 
it gives topological information that is easier to interpret and 
more directly relevant to the object of study. However, we 
need to keep in mind that with this choice, the persistence 
diagrams are not stable under the bottleneck distance. The 
results are particularly sensitive to the radius, which leaves 
holes if chosen too small and absorbs holes if chosen too big. 

The work in this paper has motivated the extension of 
Alexander duality from spaces to functions, as proved in 
0. In particular, the Euclidean Shore Theorem in that pa- 
per states that the persistence diagram of the time function 
restricted to the boundary of a Voronoi medusa can be ob- 
tained from the diagram of the time function on the Voronoi 
medusa. Generically, the boundary is a 3-manifold with- 
out boundary, so that the time function can be understood 
in Morse theoretic terms, which is sometimes convenient. 

In conclusion, we note that the framework introduced in 
this paper applies to general point processes that unfold 
in time. The latter model a variety of problems of which 



we mention a few: molecules of two fluids mixing after a 
shock-wave; microbes forming microfilms; a flock of birds 
getting into formation; two teams competing in a soccer 
game; galaxies moving under the influence of gravity. The 
measurements taken within this framework are significantly 
coarser than what could be said by studying braids [ 1 2 1 or the 
loop spaces of configuration spaces needed to detect topo- 
logical differences for particles moving in 3 dimensions JT]. 
This is not necessarily a disadvantage since the coarser in- 
formation is easier to compute as well as to comprehend. 
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